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Abstract Curved fluid interfaces are investigated on the nanometre length scale 
by molecular dynamics simulation. Thereby, droplets surrounded by a metastable 
vapour phase are stabilized in the canonical ensemble. Analogous simulations are 
conducted for cylindrical menisci separating vapour and liquid phases under con- 
finement in planar nanopores. Regarding the emergence of nanodroplets during nu- 
cleation, a non-equilibrium phenomenon, both the non-steady dynamics of conden- 
sation processes and stationary quantities related to supersaturated vapours are con- 
sidered. Results for the truncated and shifted Lennard- Jones fluid and for mixtures 
of quadrupolar fluids confirm the applicability of the capillarity approximation and 
the classical nucleation theory. 



1 Introduction 

The influence of curvature on the properties of a nanodroplet, as opposed to an 
interface that is planar (on the molecular level), is hard to capture experimentally. 
Yet it is important for refrigeration and energy technology as well as meteorology to 
understand fluid interfaces with extremely high curvatures because they characterize 
the onset of condensation and boiling processes. 

M. T. Horsch, S. K. Miroshnichenko, and J. Vrabec 

Thermodynamik und Energietechnik (ThEt), Institut fiir Verfahrenstechnik, Universitat Paderborn, 
Warburger Str. 100, 33098 Paderborn, Germany 

M. Bernreuther, C. W. Glass, and C. Niethammer 

Hbchstleistungsrechenzentrum Stuttgart (HLRS), Nobelstr. 19, 70569 Stuttgart, Germany 
M. T. Horsch, G. Jackson, and E. A. Miiller 

Molecular Systems Engineering (MSE), Centre for Process Systems Engineering, Imperial College 
London, South Kensington Campus, London SW7 2AZ, England 

* Corresponding author: Prof. Dr.-Ing. habil. Jadran Vrabec, Universitat Paderborn. E-mail: 
jadran.vrabec@uni-paderborn.de; phone: +49 5251 60 2421; fax: +49 5251 60 3522. 



1 



2 



Horsch et al. 



Beginning in the 1930s with the work of Verschaffelt [1], researchers became 
aware of the necessity of taking the internal structure of fluid interfaces into ac- 
count. They increasingly looked beyond the picture of a discrete dividing surface 
as postulated by Gibbs [2]. In the subsequent years, this led to the theoretical work 
of Guggenheim [3] and Tolman [4] which, nonetheless, was mostly still based on 
effective radii and hence on discretization. 

Today, molecular dynamics (MD) simulation provides a means of accessing the 
internal structure and the non-equilibrium behaviour of vapour-liquid interfaces di- 
rectly, on the basis of physically sound but algebraically simple effective pair po- 
tentials. For the truncated and shifted Lennard- Jones (LJ-TS) potential [5] 

Us[a l2 (r- n -r- l2 ) + a 6 (r- 6 -r- 6 )] : forr<r c , 
[ ) ~\0, for r>r c , (l) 

with a cutoff radius of r c =2.5 a, which constitutes a reliable model for the no- 
ble gases and methane [6], fluid phase boundaries have been simulated by several 
groups in the recent past [6, 7, 8, 9, 10]. Molecular simulation is particularly suit- 
able for investigating metastable states, cf. Fig. 1. Virial isotherms that accurately 
describe the conditions of fluid phase coexistence at planar or curved interfaces 

J = t bjP J \ (2) 

which were correlated to the present data, are given in Tab. 1. Therein, p is the 
pressure, T is the temperature, p is the density, and the convention k = 1 is used. 
On this basis, the present work regards both equilibrium and non-equilibrium phe- 
nomena for vapour-liquid interfaces of the LJ-TS fluid. Nucleation in supersaturated 
vapours is considered for mixtures of quadrupolar fluids as well. For a more detailed 
exposition, the reader is pointed to four recent articles [10, 1 1, 12, 13]. 



2 Curved fluid interfaces in equilibrium: Theory 

The Tolman [4] approach to curved vapour-liquid interfaces is based on comparing 
different effective radii of a droplet, namely the Gibbs adsorption radius R p , for 
which the interfacial excess density is zero, and the Laplace radius Ri from 

Y=\Rl(p,-p), 0) 

i.e. the Laplace equation in terms of the surface tension y and the pressure p, inside a 
droplet containing ; molecules under equilibrium conditions. The deviation between 
these radii, the Tolman length 

8=R P -R L , (4) 
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Fig. 1 Isotherms in a 
pressure-volume diagram 
for the LJ-TS fluid as de- 
termined by canonical MD 
simulation at temperatures of 
T = 0.65 (V), 0.75 (0), 0.85 
(A), and 0.95 e (□) as well 
as the saturated states (•) ac- 
cording to Vrabec et al. [6] in 
comparison with the present 
fifth-order virial expansion 
(— ), cf. Tab. 1 and Eq. (2). 
Results in the vicinity of the 
spinodal line are not shown 
here, and no such values were 
used for adjusting the virial 
coefficients, to ensure that 
only states unperturbed by 
nucleation were taken into 
account. 
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characterizes the curvature dependence of the surface tension [4]. The capillarity 
approximation, which postulates the droplets to be spherical and incompressible 
with a curvature independent surface tension, assumes 8 to be zero. 



Table 1 Virial coefficients for the LJ-TS fluid as determined from a fit to the MD simulation 
results shown in Fig. 1 and to the saturated vapour and liquid densities determined by Vrabec et al. 
[6]. The spinodal densities (p') and (p") for liquid and vapour, respectively, were determined 
from the virial expansion, cf. Eq. (2). 
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According to the formalism employed by Buff [14] and Kondo [15], the surface 
tension becomes mimimal if it is evaluated with respect to R^. It can be shown that 
this assertion is only valid if the interfacial area F is proportional to R^. However, 
both mechanical and thermodynamic equilibrium conditions for a droplet containing 
i molecules imply 
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) 



N,V,T 



(5) 



where V, and V are the volumes occupied by the droplet and the system as a whole, 
respectively. This only agrees with F ~ if curvature effects cancel out. 

For cylindrical interfaces, the surface tension varies with the radius according to 



an expression that is similar and analogous to Tolman's equation for droplets, where 
Ri is defined to be positive for convex and negative for concave menisci. By com- 
bining the Young equation [16] with an expansion of y to first order in 1 /R\_, the 
contact angle # of a fluid confined in a planar nanopore is obtained as 



Therein, y«, is the surface tension of the planar vapour-liquid interface, 5o» is the 
Tolman length in the planar limit, and the pore diameter is 2R m i n , while A y s indicates 
the difference between the specific surface energies of the substrate when it is in 
contact with the vapour and the liquid, respectively. 



3 Curved fluid interfaces in equilibrium: MD simulation 

Using the Isl mardyn MD program [19], equilibrium states involving droplets and 
cylindrical menisci were simulated for the LJ-TS fluid. 

Vapour-droplet equilibrium MD simulations were conducted for droplets con- 
taining on the order of 100 to 1 000 molecules in the canonical ensemble, where 
such equilibria can be stable - as opposed e.g. to the grand canonical ensemble 
where this corresponds to a free energy maximum. The droplet size was evaluated 
according to a version of the cluster criterion of ten Wolde and Frenkel [20] with the 
connectivity radius R/ : = 1.5 a and coordination numbers j > 4 defining the liquid 
phase. The present results correspond to moderately supersaturated vapours, cf. Fig. 
2, and are consistent with the results of Vrabec et al. [6] on larger droplets as well 
as the study of Napari et al. [17] covering vapours at pressures that approach the 
spinodal line. In the intermediate regime, the droplet size in equilibrium generally 
agrees well with the capillarity approximation. 

Cylindrical interfaces were investigated by simulating liquid slabs, cf. Fig. 4, 
confined between two planar and layered walls represented by coupled harmonic 
oscillators. The equilibrium positions of the wall atoms were aligned according to 
a hexagonal structure with an interatomic distance of 0.3816 (7, corresponding to 
the bond length in graphite expressed in terms of the a parameter value for fluid 
methane. Both the fluid-fluid and the fluid-wall interactions were modelled by the 




(6) 




(7) 
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LJ-TS potential, employing equal size parameters o~f w = 0~ in both cases, while the 
dispersive energy between fluid molecules and wall atoms 

£fw = (8) 

was systematically varied. The arithmetic mean of the saturated vapour and liquid 
densities was selected as a criterion for detecting the phase boundary. A circle was 
adjusted to the resulting profile at distances between 2 and 1 1 o~ from the wall, cf. 
Fig. 4, and the tangent to this circle at a distance of 1 a from the wall was examined 
to determine the contact angle. Qualitatively, the contact angles obtained by the 
present MD simulations are captured by Eq. (7), assuming a proportionality law for 

4 7s = K r (p'-p")(C-Co), (9) 

cf. Fig. 5. The magnitude of the fluid-wall dispersion for which # becomes rect- 
angular (£ = 0.1 18) was found to be temperature independent. 



4 Homogeneous vapour to liquid nucleation: Theory 

The foundations of the classical nucleation theory (CNT), concerning the first step 
of a first-order phase transition in the bulk of a metastable phase, were laid by 



Fig. 2 Droplet size ;* over 
the supersaturation ratio 
(in terms of the chemical 
potential) for vapour-droplet 
equilibria of the LJ-TS fluid 
(•) from the present work, 
(V) according to Vrabec et al. 
[6], (A) according to Napari 
et al. [17], and following the 
capillarity approximation ( — ) 
as well as the Laaksonen et al. 
[18] model (--). 
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Fig. 3 Simulation snapshots 
for the reduced fluid-wall 
dispersive energy £ of 0.09 
(left) and 0.16 (right) at a 
temperature of 0.73 e. The 
upper half is reproduced 
in the bottom to illustrate 
the effect of the periodic 
boundary condition. 





Fig. 4 Vapour-liquid interface 
profiles for the reduced fluid- 
wall dispersive energy £ of 
0.07 (A), 0.10 (□), 0.13 (o), 
and 0. 16 (V) at a temperature 
of 0.82 e. 
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Fig. 5 MD simulation results 
for r> over T with a reduced 
fluid-wall dispersive energy 
of £ = 0.07(*), 0.09 (o), 0.1 
(▼), 0.11 (□), 0.13 (■), 0.14 
(A), and 0.16 (•) as well as 
following the proportionality 
law (— , - -), cf. Eq. (7), 
with Co = 0.118, 5oo{T) from 
a correlation based on the 
data of Vrabec et al. [6], and 
a reduced fluid-wall surface 
energy difference of Ky = 7 
ere, cf. Eq. (9). 
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Volmer and Weber [21] as well as Farkas [22]. On the basis of the capillarity ap- 
proximation, the free energy of formation AA of a droplet containing ; molecules 
in the thermodynamic limit (i.e. for an infinitely large vapour phase at a constant 
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supers aturation ratio) evaluates to 

dA dF ( Pi~p\ 

^=^-(^ + Vv' (10) 

in differential terms, where fi, and p, are the chemical potential and the pressure 
inside the droplet, respectively, while jj. and p refer to the vapour. In the ; — > °° limit, 
the rightmost term of Eq. (10) yields an effective chemical potential difference 

Ati e = Li-Lk(T) + Ps( - T y P , (11) 

r 

that accounts for the 'pressure effect' [23] of the vapour - which may include the 
contribution of an inert carrier gas. Note that fi s (T) and p s (T) do not depend on 
; since these quantities characterize the saturated bulk fluid. As visualized in Fig. 
6, the presence of a carrier gas increases the free energy barrier AA* of a nucle- 
ation process, i.e. the maximum of AA reached for a critical droplet in (unstable) 
equilibrium with the vapour, corresponding to the conditions discussed above. 

From the analysis of a random walk over ;, the probability for a droplet containing 
£ molecules to eventually reach macroscopic size can be determined as 

J(^(2AA/T)dr 
Wl >~ f~exp(2AA/T)di' U } 

while the nucleation rate, i.e. the number of macroscopic liquid drops formed by 
homogeneous nucleation per volume and time unit, is 

N'zF* f-AA*\ 
j = CNT-^ r exp(^—\, (13) 



Fig. 6 Free energy of forma- 
tion according to CNT for 
CO2 droplets in supersatu- 
rated vapours at T = 250.2 K 
with S p = 2.72 and C0 2 mole 
fractions of y = 1 ( — ), 1/2 
( — ), and 1/3 (■■■), reaching 
a maximum for the critical 
droplet (□). Note that in the 
thermodynamic limit, i.e. in 
a macroscopic system, the 
Gibbs, Helmholtz, and Lan- 
dau free energies of formation 
for small droplets converge. 
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Fig. 7 Number of droplets per 
volume over simulation time 
for droplets containing ; > 10, 
25, 50, and 100 molecules 
in a canonical ensemble MD 
simulation of the LJ-TS fluid 
at T = 0.7 e and p = 0.03421 
o~~ 3 in comparison with the 
aggregated number of daemon 
interventions per volume 
in a grand canonical MD 
simulation with T = 0.7 e, 
= 2.8658, and £ = 51. 

according to CNT. In this expression, T refers to the rate at which vapour monomers 
collide with an interface (per surface area), N' is the number of monomers in the 
system, F* is the surface area of a critical droplet, and V is the system volume. The 
correction factors C, N, and z were introduced by Farkas [22], Feder et al. [24], and 
Zel'dovic [25], respectively. The overall pressure effect on J following Eq. (13), as 
discussed by Wedekind et al. [23], defines the carrier gas correction factor W. 



5 Homogeneous vapour to liquid nucleation: MD simulation 

Nucleation in supersaturated vapours was studied by simulating systems containing 
between 100 000 and 17 000 000 molecules, exploiting the excellent scalability of 
the Isl mardyn program on massively parallel computing platforms [26]. 

The method of Yasuoka and Matsumoto [27], where droplet formation rates are 
evaluated during the stage of a condensation process that corresponds to nucleation 
(rather than relaxation or droplet growth), was applied to the canonical ensemble. 
In these simulations, the vapour pressure decreased over time due to the transfer of 
molecules from the vapour to the dispersed liquid phase. Furthermore, steady state 
quantities, pertaining to nucleation only, were investigated with a new simulation 
method. This method combines the grand canonical ensemble with McDonald's 
dcemon [28], an intelligent being whose interventions eliminate droplets contain- 
ing more than I molecules; see Fig. 7 for a comparison between these approaches. 
Results for the LJ-TS fluid agree well with CNT, using a temperature independent 
value for the empirical correction factor C introduced by Farkas [22]. 

Canonical ensemble MD simulations were also conducted for multi-component 
systems containing nitrogen, oxygen, and argon - at the ratio prevalent in the earth's 
atmosphere - as well as carbon dioxide with a greater partial density than at satu- 
ration. The molecular models employed for this purpose, introduced by Vrabec et 
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Fig. 8 Nucleation rate of 
the LJ-TS fluid over super- 
saturation from the present 
MD simulations of the grand 
canonical ensemble with Mc- 
Donald's daemon (•) as well 
as according to CNT with 
C = 200 (— ) and the Laak- 
sonen et at. [18] model ( — ) 
at temperatures of T = 0.45, 
0.65, 0.7, 0.85, 0.9, and 0.95 
e. The supersaturation ratio is 
given in terms of the chemical 
potential of the vapour phase. 
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al. [29], are well-established with respect to fluid phase equilibria [29, 30, 31]. For 
these systems, the analysis of the carrier gas effect according to Wedekind et al. [23] 
is confirmed qualitatively by the determined droplet formation rates Ji, given in Tab. 
2, although significant quantitative deviations are present at high temperatures. 



6 Conclusion 

From the preceding analysis of curved vapour-liquid interfaces and homogeneous 
nucleation it can be concluded that CNT is able to capture both the nucleation rate 
and the critical droplet size for the considered systems, i.e. the LJ-TS fluid and a 
quaternary mixture of quadrupolar and unpolar fluids. The main criticism usually 
made of CNT is that it applies the capillarity approximation to small droplets where 
significant curvature effects should be expected. However, a deviation from capillar- 
ity is implicit in the prefactor C which empirically accounts for its overall influence 
on the nucleation rate. This corresponds to stating that the capillarity approximation 
overestimates the free energy barrier by T In C. 

The physical foundation of this approach is more robust than it might seem at 
first sight. By combining recent simulation results on the equilibrium vapour pres- 
sure of droplets, cf. Fig. 2, it becomes apparent that curvature effects are signi- 
ficant in the immediate vicinity of the spinodal line for the vapour, corresponding to 
;* < 100, while they are virtually undetectable for droplets containing more than 1 
000 molecules. Thus, the deviation from Eq. (10) regarding the magnitude of AA* 
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Table 2 Droplet formation rate from Yasuoka-Matsumoto (YM) canonical ensemble MD sim- 
ulation as well as critical droplet size (in molecules), Wedekind factor W, and the prediction 
J/Q CNT (£) for the droplet formation rate according to CNT with C = 1, in dependence of tem- 
perature (in units of K), supersaturation ratio (with respect to the partial density of carbon dioxide) 
and YM threshold size £ (in molecules) for the quaternary system CO2 + N2 + O2 + Ar. The mole 
fraction y of carbon dioxide in the supersaturated vapour is indicated in the table, while the com- 
position regarding the other fluids corresponds to the earth's atmosphere. The rates are given in 
units of m~ 3 s~' and where no nucleation was detected, Jem instead of //Q CNT is shown in the 
last column. 
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is dominated by an integral over the free energy of formation for extremely small 
droplets. At supersaturation ratios sufficiently distant from spinodal conditions, this 
contribution does not depend on S and can be represented by — T In C with a constant 
value of C. 
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